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Abstract 


A  new  semi-ab-initio  method  for  atomistic  simulations  based  on  the  virial  theorem  has  been 
suggested.  The  method  is  completely  within  the  realm  of  the  density-functional  theory  and  uses 
the  so-called  “reduced”  electron  spin-density  functional  (SDF).  The  crucial  component  of  the 
method  is  the  ansatz  expressing  (for  both  one-component  and  two-component  systems)  the 
electron  density  at  point  r  as  a  superposition  of  “atomic”  densities  due  to  the  neighboring  atoms. 
The  total  energy  of  this  system  is  shown  to  consist  of  three  terms.  The  first  depends  only  on  the 
simulation  volume  and  is  independent  of  the  atomic  configuration.  The  second  and  third,  like 
the  embedded-atom  method  (EAM),  are  the  interatomic  pairwise  interaction  energy,  and  an 
“N-electron”  term,  which  cannot  be  expressed  as  an  interatomic  interaction;  it  originates  from 
the  electron-correlation  interaction.  The  atomic  densities  are  constructed  using  a  set  of 
polynomial-exponential  functions  resulting  in  an  analytic  form  for  the  pair  interatomic  potential. 
The  coefficients  in  the  atomic  density  expressions  are  found  using  a  calibration  procedure  based 
on  performing  a  series  of  ab-initio  calculations  for  a  few  crystal  modifications  of  this  system. 
Success  will  depend  on  whether  the  charge  density  in  a  low-symmetry  system  under  simulation 
will  also  be  close  to  the  true  density  obtainable  from  a  meaningful  ab-initio  calculation;  then,  the 
simulation  results  would  be  identical  to  those  of  the  corresponding  ab-initio  calculation.  To 
what  extent  the  superposition  ansatz  satisfies  this  condition  is  unclear  without  experimental 
confirmation. 
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1.  Introduction 


The  formulation  of  the  density-functional  formalism  in  the  pioneering  works  of  Hohenberg 
and  Kohn  (1964),  and  Kohn  and  Sham  (1965)  published  almost  35  years  ago  opened  a  new  era  in 
the  quantum-mechanical  approach  to  the  physics  of  condensed  matter.  The  advent  of  high-speed 
computers  and  efficient  methods  of  electronic  calculations  using  the  density-functional  theory 
made  it  possible  to  achieve  tremendous  progress  in  our  understanding  of  electronic  and  structural 
properties  of  both  crystalline  and  amorphous  solids. 

The  density-functional  approach  requires  solving  a  self-consistent  Kohn-Sham 
quasi-one-particle  Schrodinger-like  equation.  While  it  is  a  rather  straightforward  and  easy  task 
for  high-symmetry  crystalline  systems,  analysis  of  low-symmetry  states,  like  crystal-lattice 
defects,  free  surfaces,  grain  boundaries  (GBs),  etc.,  is,  in  most  cases  impossible  now,  in  spite  of 
proliferation  of  high-speed  supercomputers  and  efficient  ab-initio  methods.  The  problem 
becomes  completely  intractable  when  dealing  with  arbitrary  (nonuniform)  crystal-lattice 
deformation,  defect  relaxation,  and  analysis  of  polycrystals. 

A  family  of  semi-empirical  methods*  emerged  during  the  recent  decade,  which  has  made  the 
above  problems  manageable.  In  fact,  tremendous  success  has  been  achieved  in  atomistic 
modeling  of  a  great  many  important  systems  elucidating  processes  of  atomic  relaxation 
accompanying  point  defects,  impurities,  ad-atom  layers,  GBs  and  free  surfaces,  dislocations 
dynamics,  etc.  As  a  triumph  of  atomistic  simulation,  new  and  unsuspected  features  of 
dislocation  dynamics  have  been  recently  discovered  (Bulatov  et  al.  1998). 

The  efforts  in  rectifying  semi-empirical  methods  in  recent  years  have  been  directed  at 
including  angle-dependent  potentials;  this  is  crucial  when  treating  transition  metals  and 


*  For  references  to  the  most  widely  used  method  see  the  following:  the  embedded-atom  method  (Baskes  1997);  the 
Finnis-Sinclair  N-body  potential  (Finnis  and  Sinclair  1984;  Calder  and  Bacon  1993);  the  tight-binding 
fourth-moment  method  (Carlsson  1990,  1991);  the  angle-dependent  tensor  potential,  the  so-called 
embedded-defect  methods  (Simonelly,  Pasianot,  and  Savino  1997);  and  the  bond-order  family  of  potentials 
(Tersoff  1986, 1989;  Aoki,  Horsfield,  and  Pettifor  1997;  Krasko,  Rice,  and  Yip,  to  be  published). 
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semiconductors  (Tersoff  1986,  1989;  Carlsson  1990,  1991;  Aoki,  Horsfield,  and  Pettifor  1997; 
Simonelly,  Pasianot,  and  Savino  1997).  In  attempts  to  make  semi-empirical  methods  as  close  to 
ab-initio  methods  as  possible,  the  tight-binding  approach  has  been  seriously  explored  (Bernstein 
and  Kaxiras  1996;  Yang,  Mehl,  and  Papaconstantoupolos  1998). 

Recently  (Krasko,  to  be  published),  it  was  suggested  that  the  virial  theorem,  as  applied  to  a 
condensed-matter  system,  enabled  one  to  directly  interpret  the  ingredient  contributions  to  the 
embedded-atom  method  (EAM):  the  pair  potential  and  the  embedded  function,  in  terms  of  the 
density-functional  theory.  This  interpretation  also  suggests  a  way  of  developing  a  fundamentally 
new  method  for  atomistic  simulations  that  can  be  called  “semi -ab-initio”  This  method  would  be 
using  the  virial-theorem  density-functional  expression  for  the  total  energy,  the  electron  charge 
density  being  the  only  “adjustable”  quantity.  The  method  would  be  applicable  both  to  pure 
crystals,  compounds,  and  systems  with  impurities.  Angular-dependent  interatomic  interaction 
would  also  be  automatically  included  in  a  natural  way.  The  present  technical  report  summarizes 
the  first  attempts  at  developing  such  a  method. 

The  work  summarized  in  this  technical  report  has  been  done  during  the  author’s  stay  with  the 
Department  of  Nuclear  Engineering,  MIT,  Cambridge,  MA,  as  a  visiting  scholar,  in  academic 
year  1997-1998. 

2.  Theoretical  Background  of  the  New  Method 

As  mentioned  above,  in  its  traditional  implementation,  the  density-functional  approach 
requires  solving  a  self-consistent  Kohn-Sham  quasi  one-particle  Schrodinger  equation.  As  for 
the  exact  density  functional,  or  the  spin-density  functional  (SDF)  for  the  total  energy,  E{p(r)} 
[where  p(r)  is  the  electron  density  at  point  r],  its  explicit  form  is  not  known  and  may  never  be. 

In  the  recent  years,  the  interest  toward  finding  a  meaningful  approximation  for  the  SDF, 
using  various  versions  of  “orbital  free”  kinetic  energy  (KE)  functionals  has  been  awakened  (e.g., 
Wang  and  Teter  1992;  Pearson,  Smargiassi,  and  Madden  1993;  Smargiassi  and  Madden  1994). 
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On  the  other  hand,  the  SDF,  in  what  can  be  called  a  “reduced”  form,  can  easily  be  found 
using  the  scaling  procedure  usually  involved  in  formulating  the  virial  theorem  (Ross  1969).  In 
this  section,  we  show  how  such  a  reduced  SDF  can  be  found. 

The  total  energy  functional,  E{p(r)>,  is,  as  usual, 

E{p(r)}  =  T{p(r)}  +  U{p(r)>  +Exc{p(r)> ,  (1) 

where  T{p(r)>,  U{p(r)},  and  ExC{p(r)}  are,  respectively,  the  KE  of  noninteracting  electrons  and 
the  potential  and  exchange-correlation  energy  functionals.  Note  that  here,  p(r)  may  be  a  trial 
density. 

The  potential  energy  functional  is 

U  =  J  V(r)  p(r)  dr  +  1/2  J  J  p(r)  p(r')/lr-r'l  dr  dr'  +  Ees ,  (2) 

where  V(r)  is  the  “external”  (electron-ion)  potential,  and  Ees  is  the  electrostatic  ion-ion 
interaction  energy.  As  for  ExC{p(r)},  we  suppose  that  it  is  known  either  in  a  local  density 
approximation  (LDA)  (Hohenberg  and  Kohn  1964)  or  the  so-called  “generalized  gradient 
approximation”  (GGA)  (Perdew,  Burke,  and  Emzerhof  1996): 

Exc  =  J  p(r)  £Xc(r)  dr ,  (3) 

where  exc(r)  =  £xc(p(r))  in  LDA  and  £xc(p(r),  t,  s )  in  GGA  (where  t  and  s  are  the  dimensionless 
density  gradient  variables). 

Following  the  scaling  procedure  used  by  Ross  (1969),  let  us  introduce  the  dimensionless 
coordinates  x  =  r/Q1/3. 
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Then, 


P(r)  ->  p(x)/Q , 

U(p(r))->  U{p(x)}/£21/3 ,  and 

ExC{p(r) >  -»  Q  J  (P(t)/Q)  exc(T)  dx  (4) 

(the  GGA  dimensionless  gradient  variables  t  and  s  are  independent  of  volume).  The  KE  of 
noninteracting  electrons  is 

T(x)  =  ZoccJ  (^(x)*  (1/2)  V2<tH(x)  dx ,  (5) 

where  <(h(x)’s  are  the  (trial)  one-electron  wave  functions  and  the  summation  is  extended  to  all 
occupied  states.  Hence,  T  scales  as 

T{p(r)>  ->  T(x)/QM.  (6) 

Now,  let  us  calculate  P  =  -dE  {p(r)}/d£2: 

P  =  2/3  T/Q  -  d  (U  +  Exc)/dQ  .  (7) 

The  volume  derivative  is  easily  taken  using  the  scaling  equations  (4)  and  (6): 
d(U  +  Exc)/dQ  =  -U(Qy(3Q) 

+1  /Q  J[p(r)  exc(r)  -  p+(r)  M*c+(r)  -  p'  (r)  jixC"  (r)]dr  (8) 
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If  p(r)  is  the  exact  solution  of  the  Kohn-Sham  equation  (the  minimizer  of  the  SDF),  then 
P  =  P,  the  physical  pressure,  and  equation  (7)  is  just  one  of  the  versions  of  the  virial  theorem.  It 
may  be  represented  in  the  traditional  form  (Slater  1972;  Janak  1974)  using  equation  (8): 

2T  =  3  PQ  -  U  +  3  J  [p(r)  exc(r)  -  p+(r)  p*c+(r)  -  p"  (r)  p*c'  (r)]dr .  (9) 

In  equation  (8)  and  here, 

|ixc+,~  (r)  =  8Exc/8p+’“(r)  (10) 

is  the  exchange-correlation  one-electron  potential  Cr  stands  for  spin  polarization).  The 
expression  for  the  exchange-correlation  part  of  equation  (9),  here,  is  more  general  than  that  for 
non-spin-polarized  systems  in  Janak  (1974).  Since  the  GGA  gradient  variables  do  not  depend  on 
volume,  equation  (9)  is  also  the  expression  for  the  virial  theorem  in  GGA. 

We  again  stress  that  both  the  scaling,  equations  (4)  and  (6),  and  the  taking  of  the  derivative 
over  the  volume  are  valid  not  for  just  the  SDF  of  stationary  electron  density,  p(r)  but  for  any 
appropriate  trial  density.  In  what  follows,  p(r)  means  just  such,  a  trial  density. 

It  should  be  also  stressed  that,  although  the  virial  theorem,  equation  (9),  is  also  valid  for  trial 
densities  (then  P  has  to  be  substituted  for  P),  it  is  satisfied  for  a  real  physical  pressure  P  only  if 
p(r)  is  a  highly  accurate  solution  of  the  corresponding  Kohn-Sham  equation. 

Let  us  return  to  equation  (7).  We  have 

T  =  3/2PQ  +  3/2Qd(U  +  Exc)/d£2.  (11) 


If  P  =  0,  then 


T  =  3/2  Q0  d(U  +  ExcVdQlo^o , 


(12) 


where  Qo  is  the  volume  corresponding  to  P  =  0. 

For  an  arbitrary  volume,  £2,  the  noninteracting  electron  KE,  due  to  the  scaling,  equation  (6), 
equals 

T(Q)  =  [3/2  Q0  d(U  +  ExC)/dQb=fx>  ]  (Qq/Q f3  (13a) 

or,  making  use  of  equation  (8), 

T(Q)  =  WQ)m{-1/2  U(Q0) 

+  3/2  J[p(r)  exc(r)  -  p+(r)  M*c+(r)  -  p  "  (r)  p*c '  (r)]drb=o0> .  (13b) 

Finally,  the  total  energy  functional  is 

E(Q)  =  U  +  Exc  +  3/2  Q  (Qo/Qf3  d(U  +  ExC)/dQb=Qo  (14a) 

or,  in  the  explicit  form, 

E(Q)  =  U  +  Exc  +  (Qo /Qf3{- 1/2  U(Qo)  +  AJoo,}  ,  (14b) 

where 


Axc  =  3/2  J[p(r)  £xc(r)  -  p+(r)  Pxc+(r)  -  p'  (r)  M*c‘  (r)]dr .  (15) 

Then,  the  “pressure”  equals 

P(Q)  =  -d(U  +  Exc)/dQ  +  (Qo/Q)5/3d(U  +  ExC)/dQb=Oo  (16a) 
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or 


P(Q)  =  1/(3Q){-U(Q)  +  2  Axc  -  (Qo/nf3  [-U(Q0)  +  2  Axclo=n0]  >  •  (16b) 

Equations  (14a)  and  (14b)  are  the  reduced  SDF  that  we  sought.  It  is  an  explicit  functional  of 
p(r)  valid  for  any  appropriate  trial  density.  For  a  stationary  density  (satisfying  a  Kohn-Sham 
equation),  equations  (7)  and  (9)  are  equivalent  to  the  virial  theorem  and  equations  (14a)  and 
(14b)  just  “recalculate”  the  total  energy,  which  has  to  be  equal  to  that  following  from  the 
traditional  Kohn-Sham  method.  Equation  (16)  is  then  the  equation  of  state. 

Of  course,  the  reduced  SDF  does  not  contain  all  the  information  of  the  general  universal 
SDF.  However,  it  may  serve  as  a  basis  for  a  fundamentally  new  approach  toward  developing 
s&mx-ab-initio  methods  of  atomistic  simulations. 

In  the  next  section,  we  outline  the  algorithm  of  constructing  the  model  density  and  the 
procedure  of  calculating  the  total  energy,  which  lies  at  the  core  of  the  new  proposed  method. 

3.  Formulation  of  the  Method 

Suppose  that  the  ground-state  electron  density  may  be  written  as  a  superposition  of  “atomic” 
(or  “pseudo-”  or  “quasi-atomic”)  densities:* 

p(r)  =  XRp°(lr-RI).  (17) 

This  ansatz  has,  in  the  past,  been  studied  in  great  detail.  It  was,  in  particular,  shown  (Chetty, 
Jacobsen,  and  Norskov  1991)  that  optimized  and  transferable  atomic  densities  can  be  found  from 


*  There  is  no  unique  and  exact  way  of  subdividing  the  total  electron  charge  density  into  a  superposition  of  atomic 
densities.  In  fact,  by  introducing  Wannier  functions,  one  can  show  that,  in  general,  die  electron  charge  density 
cannot  be  represented  exacdy  as  the  superposition  of  neighbor-centered  positive  functions  that  could  be  interpreted 
as  individual  atomic  densities;  it  is  possible  only  for  insulators  (Krasko,  to  be  published). 
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the  first  principles.  Equation  (17)  is  the  main  approximation  of  our  method.  A  procedure  of 
constructing  p°(lr  -  Ri)  is  discussed  next.  The  ansatz,  equation  (17),  also  enables  one  (Krasko 
1999)  to  interpret  the  potential  energy  contributions  to  E°,  equation  (10),  in  an  intuitively 
appealing  form,  the  one  similar  to  that  of  the  EAM. 

Using  equation  (17),  the  potential  energy,  U,  equation  (2),  can  be  written  down  in  terms  of  an 
effective  pair  potential: 

1/2U  =  E°  +  1/2  2W  F(R  -  R') .  (18) 

The  first  term  depends  only  on  the  volume  but  not  the  atomic  coordinates  (we  give  its  explicit 
expression  later).  The  second  term  is  the  pair  interaction  energy  with  the  pair  potential: 

V(R)  =  - 1/(271?)  Jdq  Z  p°(q)/q2  exp(iq-R) 

+  1/(471?)  J dq  lp°(q)l2/q2  exp(iq-R)  +  Z2/IRI .  (19) 

Here  and  below,  p°(q)  is  the  Fourier  transform  of  the  atomic  charge  density,  p°(r)  and  the 
integration  extends  over  all  the  reciprocal  space. 

The  exchange-correlation  contribution,  Axc,  to  the  total  energy,  ExC,  from  equation  (3), 

Axc  =  3/2 J  p(r)  [e(r)  -  p(r)]  dr,  (20) 

cannot,  however,  be  represented  in  terms  of  an  interatomic  potential. 

Finally  the  total  energy,  equation  (14b),  can  be  rewritten  as 

E  =  E°+  1/2  (1  -  1/2  (£VG)1/3)W  V(R  -  R') 

+  Exc  +  (Qo/af3  Axc{p(r)  >h=no .  (21) 
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Here, 


E°  =  -N  (1  -  1/2  (Qo/Q)1/3)  {1/(4 7C2)  Z  Jdq  p°(q)/q2 

+  l/CSTt2)  Jdq  lp°(q)l2/q2>  (22) 

(N  is  the  total  number  of  atoms  in  the  simulation  volume).  The  origin  of  this  energy  contribution 
is  obvious;  it  is  the  atomic  “self-interaction,”  the  term  corresponding  to  R  =  R\  which  has  been 
excluded  from  the  pair  potential  sum  in  equation  (20). 

Now,  we  proceed  with  formulation  of  the  suggested  method.  As  was  mentioned  previously, 
our  goal  is  to  parametrize  the  atomic  density,  p°(r). 

Since  it  is  desirable  to  represent  the  density  as  exactly  as  possible,  it  is  appealing  to  separate 
the  contributions  of  the  core  and  the  valence  electrons: 

p°(lrl)  =  p°cor(lrl)  +  p°vai(lrl) .  (23) 

p°cor(lrl)  will  be  approximated  once  and  for  all  for  the  given  atom  using  the  results  of  ab-initio 
calculations.  Let 


P°cor(lrl)  =  AC0rEI-{c°i+  c\r  +  c2  r2  +...)  exp(-X.c;r)>.  (24) 

Here,  Ac0r  is  the  normalization  constant  (securing  the  right  number  of  core  electrons),  and  cVs 
and  A. Vs  are  the  parameters  to  be  found  by  root-mean-square  (rms)-fitting  of  equation  (24)  to  the 
calculated  core  densities  (depending  on  the  core  structure,  there  may  be  a  few  different 
“subshells,”  i). 

The  valence  charge  density  will  be  approximated  by  the  same  class  of  functions: 


P°vai(lrl)  =  AvaiEi{v°,+  v^  r  +  v2  r2  +...)  exp(-yv,  r)> , 


(25) 
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where  Avai,  again,  is  the  normalization  constant  (securing  the  right  number  of  valence  electrons) 
and  vY  s  and  yY  s  are  the  adjustable  parameters.  In  a  transition  metal,  one  may  expect  three 
different  subshells,  i,  corresponding  to  -s,  -p,  and  -d  electrons. 

The  chosen  class  of  functions  is  convenient,  since  it  allows  the  integration  to  be  performed 
analytically  when  calculating  both  the  Fourier  transforms  of  the  densities  and  the  effective  pair 
potential,  equation  (19),  and  the  structure-independent  term,  E°,  equation  (22).  (We  are  not 
writing  down  those  lengthy  formulae  here). 

Thus,  the  first  two  contributions  to  E,  equation  (21),  can  easily  be  found  in  analytic  form. 
The  calculation  of  the  third  and  fourth  terms  to  do  with  exchange  correlation  is  less 
straightforward.  First  of  all,  an  exchange-correlation  model  (either  local  or  GGA)  has  to  be 
chosen  among  a  few  models  being  currently  used  in  ab-initio  methods.  The  corresponding 
subroutine  [having  p(r)  as  an  input  and  e(r)  and  |i(r)  as  outputs]  may  be  borrowed  from  any 
ab-initio  code  in  the  public  domain.  What  is  important  is  that  the  same  exchange-correlation 
approximation  be  used  both  in  evaluating  Axc  and  all  the  ab-initio  calculations  used  in  the 
method  calibration.  We  suggest  to  use  the  most  advanced  exchange-correlation  model,  GGA 
(Perdew,  Burke,  and  Emzerhof  1996). 

The  calculation  of  the  exchange-correlation  contribution  in  equation  (21),  should  be 
performed  as  follows.  Within  the  “simulation  volume,”  a  grid  of  coordinates  {r}  is  to  be 
constructed.  Based  on  our  experience  with  ab-initio  calculations,  a  nonuniform  grid  has  to  be 
used,  with  the  grid  points  thickening  around  the  lattice  cites,  {R>,  and  being  looser  in  the 
interstitial  volume.*  Then,  the  total  density,  p(r),  as  a  superposition  of  p°(lr  -  Rl)’s  is  calculated 
on  this  grid.  Computationally,  this  is  to  be  a  rather  fast  procedure  since,  typically,  p°(lr  -  Rl) 
virtually  completely  decays  by  the  third  shell  of  neighbors.  As  the  last  step  in  the  calculation  of 
the  exchange-correlation  contributions  to  the  total  energy,  the  integrations  are  to  be  performed 
using  one  of  the  available  methods  of  numerical  three-dimensional  (3-D)  integration. 


*  Recently,  a  more  sophisticated  way  of  introducing  an  “adaptive”  coordinate  system  was  suggested;  see  Modine, 
Zumbach,  and  Kaxiras  (1997). 
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The  main  idea  behind  the  reduced  SDF  is  to  scale  the  KE  functional  and  thus  make  use  of  the 
virial  theorem.  As  a  result,  the  knowledge  of  the  electron  charge  density  (as  well  as  both  the 
potential  and  the  exchange  correlation  contributions)  at  P  =  0,  or  P  =  0  becomes  a  necessity. 
Therefore,  in  every  calculation  of  the  total  energy  at  a  given  atomic  configuration  and  simulation 
volume,  a  part  of  the  calculation  must  be  performed  at  an  “equilibrium”  volume,  Q0.  This 
volume,  however,  is  unknown. 

We  suggest  the  following  procedure  of  finding  Q<>  The  total  energy  at  volume  Q, 
equation  (21),  depends  on  Q0  •  We  may  consider  Q0  as  a  variational  parameter  and  minimize 
equation  (21)  at  fixed  Q,  with  respect  to  Q0  •  The  obtained  value  will  also  automatically  satisfy 
equation  (16b)  at  P  =  0. 

In  order  to  deal  with  both  an  arbitrary  volume  and  Qo,  we  again,  use  a  scaling  procedure.  Let 
Qr  be  a  “reference”  simulation  volume.  Actually,  Qr  can  be  any  volume  within  the  range  of 
physically  meaningful  volumes.  Then  the  arbitrary  volume  can  be  represented  as 


Q  =  C53Qt  (26) 

where  03  is  a  scaling  factor.  Then,  all  the  distances  in  our  model  will  be  scaled  uniformly,  and 
the  density,  equation  (17),  corresponding  to  the  volume,  Q,  is  found  simply  by  multiplying  all  the 
r’s  and  R’s  by  gj.  Suppose  Q  is  the  volume  at  which  the  total  energy,  equation  (21),  is  to  be 
calculated.  The  equilibrium  volume,  Qo,  is  then  the  volume  corresponding  to  the  minimum  of 
the  total  energy,  equation  (21),  at  fixed  Q.  Obviously,  it  would  correspond  to  the  scaling  factor 
G5o: 


Q0  =  G5o3&  .  (27) 

Thus,  we  have  found  the  equilibrium  volume,  Qo-  From  this  point  on,  the  procedure  of 
calculating  the  total  energy,  equation  (21),  is  straightforward. 
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This  completes  the  procedure  of  calculating  the  total  energy  of  the  system. 

We  now  summarize  the  method  outlined  previously.  Suppose  we  want  to  calibrate  the 

method  for  transition  metal  Mo.  The  following  steps  will  have  to  be  followed. 

(1)  We  begin  by  performing  a  series  of  ab-initio  calculations  for  body-centered  cubic  (bcc)  (e.g., 
FLAPW  [Singh  1994])  to  find  Ebcc(Q)  and  Q0bcc  (P  =  0).  Repeat  the  previous  calculations 
for  face-centered  cubic  (fee),  hexagonal  close-packed  (hep),  and,  possibly,  a  low-symmetry 
hypothetical  phase  (HP)  of  Mo,  and  find  Efcc(Q),  Ehcp(Q),  EHP(Q)  and  the  P  =  0  volumes: 
Q0icc,  Q0hcp,  Oo®.  These  calculations  also  generate  the  core  charge  densities  to  be  used  in 
the  next  step. 

(2)  Approximate  the  Mo  core  electron  density  using  equation  (24)  and  find  the  appropriate 
parameters  by  rms-fitting  the  function  to  the  calculated  ab-initio  Mo  core  density. 

(3)  After  the  previous  ab-initio  calculations  have  been  performed,  we  have  six  equations  for 
fitting  the  parameters  of  the  valence  electron  density,  equation  (24):  three  for  the  P  =  0 
energies  of  the  three  phases  [we  need  to  calculate  only  E,  equation  (21)]  and  three  equations 
P(Q0)  =  0.  Some  more  equations  may  be  used;  for  example,  fitting  the  total  energy  to  the 
ab-initio  energies  for  other  HPs  or  the  elastic  moduli:  Cn,  C12,  C44,  and  B  (we  do  not  write 
out  the  corresponding  formulae,  as  well  as  the  expression  for  the  force  on  each  atom;  they 
are  quite  straightforward  but  require  a  tedious  algebra). 

(4)  Having  solved  this  system  of  simultaneous  equations,  one  will  find  the  coefficients  in  the 
approximation,  equation  (25).  The  more  complicated  procedure  of  finding  £2o  will  be 
necessary  for  atomistic  simulations  outside  the  ab-initio  calibration. 

(5)  Now,  with  all  the  coefficients  in  equations  (24)  and  (25)  known,  the  calibration  of  the 
method  is  completed  and  one  may  proceed  with  atomistic  simulations  of  the  system  of 
interest. 


f 
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As  an  example  of  such  a  system,  let  us  consider  a  grain  boundary.  Let  the  simulation  volume 
(for  the  unrelaxed  GB)  be  Q.  that  contains  N  atoms.  The  Mo  atoms  take  positions  {R>.  First,  we 
calculate  £20  by  minimizing  E,  equation  (21),  at  the  fixed  simulation  volume  Q,  with  respect 
to  Q0- 


In  order  to  analyze  the  atomic  relaxation,  one  has  to  calculate  the  forces  on  each  atom,  then 
vary  the  atomic  positions  {R}  toward  the  state  of  mechanical  equilibrium  (forces  equal  zero),  and 
repeat  calculations  until  an  absolute  energy  minimum  has  been  achieved.  On  each  step  of 
shifting  the  atomic  positions,  the  new  O0  has  to  be  calculated  using  the  previously  outlined 
procedure. 

Here,  we  should  again  stress  that  the  total  energy,  equation  (21),  does  have  the  variational 
properties  with  respect  to  the  electron  charge  density.  Therefore,  the  minimization  to  be 
performed  does  mean  finding  the  density  that  would  correspond  to  the  energy  minimum,  albeit  in 
a  restricted  manifold  of  the  superpositions  of  “atomic  densities,”  equation  (17).  Thus,  the 
precision  of  the  method  as  a  whole  depends  chiefly  on  the  generation  of  realistic  density, 
equation  (17). 

In  the  next  section,  we  show  how  the  method,  in  a  natural  way,  can  be  generalized  to  a  two- 
component  system:  a  crystal  with  impurity  atoms  or  a  compound. 

4.  Generalization  of  the  Method  to  Two-Component 
Systems 

The  general  formalism  of  the  virial  theorem  is  valid,  irrespective  of  the  system  it  is  applied 
to.  The  difference  between  a  one-component  and  a  two-component  (or,  for  that  matter,  a 
multicomponent)  system  is  only  in  how  the  electron  density,  p(r),  is  represented. 
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Let  the  system  be  a  substitutional  solid  solution  of  components  A  and  B  (the  generalization  to 
interstitial  solid  solution  is  straightforward).  Then,  the  electron  charge  density  at  a  given  point 
can  be  written  as  a  generalized  ansatz: 

p(r)  =  X  R  {CA(R)p0A(r  -  R)  +  CB(R)p°B(r  -  R» .  (28) 

The  “external,”  ion-electron  potential  is 

V(r)  =  X  r{Ca(R)  Za  +  Cb(R)  Zb  }/lr  -  Rl,  (29) 

where  CA(R)  and  CB(R)  are  the  random  quantities  taking  the  value  1  if,  respectively,  an  atom  A 
or  an  atom  B  sits  in  the  site  R;  and  0  otherwise  (CA(R)  +  CB(R)  =  1),  p°A(r  -  R),  and  p°B(r  -  R) 
are,  as  before,  the  atomic  densities  and  ZA  and  ZB  are  the  nuclei  charges  of  atoms  A  and  B. 

Then  E  [a  counterpart  of  equation  (21)]  can  be  written  in  terms  of  effective  pair  potentials  as 

E  =  E°+  1/2  (1  -  1/2  (Qo/Q)1/3)Xr*r.  {Ca(R)Ca(R')Vaa(R  -  R') 

+  CB(R)CB(R')FBB(R  -  R')  +  2  Ca(R)Cb(R')Vab(R  -  R')> 

+  E*c  +  (Q0/Q)2/3Axc{p(r)  >lfi=oo ,  (30) 

where  ExC  and  Axc{p(r)>  are  defined  by  the  same  expressions,  equations  (3)  and  (15)  but  with  the 
density  p(r)  obeying  the  new  ansatz,  equation  (28). 

The  first  term,  again,  depends  only  on  volume: 

E°  =  -  l/(47t2)  <Z>Jdq  <p°(q)>/q2  +  1/(87^)  Jdq  <lp°(q)l2>/q2  ,  (31) 


where 


<Z>  =  cAZA  +  cBZB, 
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<P°(q)>  =  cAp°A(q)  +  cBp°B  (q) , 


and 


<lp°(q)l2>  =  cAlp°A(q)l2  +  cBlp°B(q)l2  •  (32) 

Also,  cA  and  Cb  are  the  atomic  fractions  of  atoms  A  and  B: 


ca=SrCa(R)/N 


and 


Cb=  ZrCb(R)/N, 


where  N,  again,  is  the  total  number  of  atoms  in  the  simulation  volume.  The  effective  interatomic 
potentials  then  are  defined  as  (a,  (5  =  A,  B): 

^ap(R)  =  -Za/(27t2)  J dq  p°p(q)/q2  exp(iq-R) 

+  l/(4it2)  | dq  p°a(q)  p°p(q)*/q2  exp(iq-R)  +  Za  Zp/IRI .  (33) 

The  formulation  of  the  method  for  two-component  systems  is  virtually  the  same  as  that 
outlined  before.  Now  one  has  two  core  densities:  p°’A’BCor,  and  two  valence  densities:  p0’A,Bvai- 
The  calibration  and  implementation  of  the  method  is  similar  to  those  in  one-component  case;  the 
only  difference  being  that  the  number  of  fitting  equations  doubles  and  the  equations  have  to 
reflect  some  representative  high-  and  low-symmetry  two-component  configurations. 

5.  Conclusion 

To  summarize,  we  suggest  a  semi  -ab-initio  method  for  atomistic  simulations  based  on  the 
virial  theorem.  The  method  is  completely  within  the  realm  of  the  density-functional  theory.  The 
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crucial  component  of  the  method  is  ansatz  equation  (17)  [or  equation  (28)  for  a  two-component 
system]  expressing  the  electron  density  at  point  r  as  a  superposition  of  atomic  densities  due  to 
the  neighboring  atoms.  If,  upon  calibrating  the  method,  the  E’s  are  fitted  to  the  equilibrium 
energies  of  a  few  crystal  modifications,  the  electron  densities  found  by  the  atomic  densities 
superposition  for  those  modification  are  very  close  to  the  “true”  stationary  densities  that  could  be 
found  from  an  ab-initio  calculation.  Thus,  the  success  of  the  method  will  depend  on  whether  the 
charge  density  in  a  low-symmetry  system  under  simulation  is  also  close  to  the  true  density  that 
could  be  obtained  from  a  meaningful  ab-initio  calculation;  then,  the  results  of  the  simulation 
would  be  identical  to  those  of  the  corresponding  ab-initio  calculation.  To  what  extent  the 
superposition  ansatz  satisfies  this  condition  is  unclear,  a  priori  testing  and  experimentation  are 
necessary.  As  was  mentioned  previously,  according  to  Chetty,  Jacobsen,  and  N0rskov  (1991), 
meaningful  transferable  ab-initio  atomic  densities  could  be  found. 

Regarding  the  separation  of  core  and  valence  densities  [equations  (23)-(25)],  a  note  should 
be  added.  “Fixing”  the  core  density  is  reminiscent  of  a  so-called  “frozen  core”  approximation. 
When  the  frozen-core  approximation  is  involved  in  an  ab-initio  method,  the  virial  theorem  is  not 
valid.  The  reason  for  this  failure  is  that,  while  the  frozen  core  potential  is  treated  as  a  part  of  the 
external  potential,  the  exchange-correlation  core-valence  interaction  is  also  preserved  (the  virial 
theorem  would  be  valid  again  if  the  latter  were  dropped).  In  our  case,  however,  the  frozen  core  is 
a  part  of  the  total  electron  density  and  it  is  the  total  electron  density  that  is  subject  to  “fitting”  to 
the  ab-initio  density.  In  fact,  from  the  physical  point  of  view,  it  would  probably  be  more  logical 
if  the  calibration  of  the  method  consisted  in  directly  fitting  the  model  density  to  the  ab-initio 
density  rather  than  requiring  that  the  calculated  and  ab-initio  energies  were  equal.  We  are  going 
to  explore  this  path  also  in  a  subsequent  paper. 

As  was  already  mentioned,  the  new  method  is  quite  versatile.  If  additional  ab-initio 
information  about  the  system  of  interest  is  available,  it  can  be  easily  implemented  to  make  the 
method  more  realistic.  As  an  example,  while  doing  simulations  on  diamond-type 
semiconducting  crystals,  it  would  be  helpful  (perhaps  even  necessary)  to  introduce,  as  a  separate 
electron  density,  the  so-called  “bond  charges”  (Phillips  1973).  The  method  can  also  be  applied  to 
atomistic  simulation  of  ferromagnetic  systems;  in  that  case,  the  separate  spin-up  and  spin-down 
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densities  have  to  be  introduced  (the  exchange-correlation  contributions  to  the  energy  can  be 
easily  generalized  for  a  spin-polarized  case). 

An  important  feature  of  the  method  is  that  it  automatically  allows  for  “directional  bonds.”  In 
spite  of  the  fact  that,  in  the  suggested  implementation,  the  total  valence  charge  is  a  superposition 
of  spherically  symmetric  charges,  p(r)  does  have  the  symmetry  of  the  crystal  lattice,  thus 
incorporating  all  “directionality”  features.  It  is  interesting,  however,  that  the  effective  pah- 
potential,  equation  (19)  (and  its  counterparts  for  two-component  system),  is  lacking 
directionality.  Since  the  Fourier  transform  of  the  density,  p(q)  =  p(lql),  the  angle  integrations  in 
equation  (19)  are  straightforward,  resulting  in  the  remaining  integral  in  equation  (19)  having  the 
integrand  that  is  proportional  to  sin(qR)/(qR),  where  q  and  R  are  the  moduli  of  the  corresponding 
vectors.  As  a  result,  the  pah  potential  depends  only  on  the  distance: 

V(R)  =  V(IRI) . 

hi  fact,  that  is  what  is  to  be  expected  of  a  pah  potential;  if  the  space  is  isotropic,  the  two 
atoms  interacting  via  a  pah  potential  do  not  have  to  know  anything  about  their  environment  since 
no  special  dhection  has  been  singled  out. 

The  directionality,  however,  is  completely  preserved  in  the  exchange-correlation 
contribution,  Axc{p(r)>,  equation  (15),  through  the  symmetry  of  p(r). 

As  was  previously  mentioned,  the  method  is  quite  versatile,  and,  if  necessary,  nonspherically 
symmetric  atomic  densities  may  be  introduced  (Lowdin  1956). 


17 


Intentionally  left  blank. 


18 


6.  References 


Aoki,  M. ,  A.  Hdrsfield,  and  D.  G.  Pettifor.  “Tight-Binding  Bond  Order  Potential  and  Forces  for 
Atomistic  Simulations.”  Journal  of  Phase  Equilibria,  vol.  18,  p.  614, 1997. 

Baskes,  M.  I.  “Determination  of  Modified  Embedded  Atom  Method  Parameters  for  Nickel.” 
Materials  Chemistry  and  Physics,  vol.  50,  p.  152, 1997. 

Bernstein,  N.,  and  E.  Kaxiras.  “Nonorthoganal  Tight-Binding  Hamiltonian  for  Defects  and 
Interfaces  in  Silicon.”  Physical  Review,  vol.  B  56,  p.  10488, 1996. 

Bulatov,  V.,  F.  Abraham,  L.  Kubin,  B.  Devincre,  and  S.  Yip.  “Dislocation  Junctions  and  Crystal 
Plasticity:  Linking  Atomistic  and  Mesoscale  Simulations.”  Nature,  vol.  391,  p.  669, 1998. 

Calder,  A.  F.,  and  D.  J.  Bacon.  “A  Molecular  Dynamics  Study  of  Displacement  Cascades  in 
a-Iron.”  Journal  of  Nuclear  Materials,  vol.  207,  p.  25,  1993. 

Carlsson,  A.  E.  “Beyond  Pair  Potentials  in  Elemental  Transition  Metals  and  Semiconductors.” 
Solid  State  Physics,  vol.  43,  p.  1,  H.  Ehrenreich  and  D.  Turnbull  (editors),  1990. 

Carlsson,  A.  E.  “Angular  Forces  in  Group- VI  Transition  Metals:  Application  to  W(100).” 
Physical  Review,  vol.  B  44,  p.  6590, 1991. 

Chetty,  N.,  K.  W.  Jacobsen,  and  J.  Norskov.  “Optimized  and  Transferable  Densities  From 
First-Principles  Local  Density  Calculations.”  Journal  of  Physics:  Condensed  Matter,  vol.  3, 
p.  5437, 1991. 

Finnis,  M.  W.,  and  J.  E.  Sinclair.  “A  Simple  Empirical  N-body  Potential  for  Transition  Metals.” 
Philosophical  Magazine,  vol.  A50,  p.  45, 1984  (errata,  vol.  A53,  p.  161, 1986). 

Hohenberg,  P.,  and  W.  Kohn.  “Inhomogeneous  Electron  Gas.”  Physical  Review,  vol.  136, 
p.  864, 1964. 

Janak,  J.  F.  “Simplification  of  Total-Energy  and  Pressure  Calculations  in  Solids.”  Physical 
Review,  vol.  B  9,  p.  3985, 1974. 

Kohn,  W.,  and  L.  R.  Sham.  “Self-Consistent  Equations  Including  Exchange  and  Correlation 
Effects.”  Physical  Review,  vol.  A  140,  p.  1133, 1965. 

Krasko,  G.  L.  “The  Virial  Theorem  and  the  Embedded  Atom  Method.”  To  be  published. 

Krasko,  G.  L.,  B.  Rice,  and  S.  Yip.  “A  Bond-Order  Potential  for  Atomistic  Simulations  in  Iron.” 
To  be  published. 


19 


Lowdin,  P.  O.  “Quantum  Theory  of  Cohesive  Properties  in  Solids.”  Advances  in  Physics, 
vol.  64,  chaps.  2.1.3  and  6.1,  p.  1, 1956. 

Modine,  N.  A.,  G.  Zumbach,  and  E.  Kaxiras.  “Adaptive-Coordinate  Real-Space  Electronic 
Structure  Calculations  For  Atoms,  Molecules,  and  Solids.”  Physical  Review,  vol.  B  55, 
p.  10289, 1997. 

Pearson,  M.,  E.  Smargiassi,  and  P.  A.  Madden.  “Ab  initio  Molecular  Dynamics  With  an 
Orbital-Free  Density  Functional.”  Journal  of  Physics:  Condensed  Matter,  vol.  5,  p.  3221, 
1993. 

Perdew,  J.  P.,  K.  Burke,  and  M.  Emzerhof.  “Generalized  Gradient  Approximation  Made 
Simple.”  Physical  Review  Letters,  vol.  77,  p.  3865, 1996. 

Phillips,  J.  C.  Bands  and  Bonds  in  Semiconductors.  Academic  Press,  NY,  1973. 

Ross,  M.  “Pressure  Calculations  and  the  Virial  Theorem  for  Modified  Hartree-Fock  Solids  and 
Atoms.”  Physical  Review,  v ol.  179,  p.  612, 1969. 

Simonelly,  R.,  R.  Pasianot,  and  E.  J.  Savino.  “Phonon  Dispersion  Curves  for  Transition  Metals 
Within  the  Embedded-Atom  and  Embedded-Defect  Methods.”  Physical  Review,  vol.  B  55, 
p.  5570, 1997. 

Singh,  D.  J.  Plane  Waves,  Pseudopotentials  and  the  LAPW  Method.  Kluger  Academic 
Publishing,  Boston,  MA,  1994 

Slater,  C.  “Hellman-Feynman  and  Virial  Theorems  in  the  X-a  Method.”  Journal  of  Chemical 
Physics,  vol.  57,  p.  2389, 1972. 

Smargiassi,  E.,  and  P.  A.  Madden.  “Orbital-Free  Kinetic  Energy  Functionals  for  First-Principles 
Molecular  Dynamics.”  Physical  Review,  vol.  B  49,  p.  5220, 1994. 

Tersoff,  J.  “New  Empirical  Model  for  Structural  Properties  of  Silicon.”  Physical  Review 
Letters,  vol.  56,  p.  632, 1986 

Tersoff,  J.  “Modeling  Solid-State  Chemistry:  Interatomic  Potentials  for  Multicomponent 
System.”  Physical  Review  ,  vol.  B  39,  p.  5566, 1989. 

Wang,  L.-W.,  and  M.  P.  Teter.  “Kinetic-Energy  Functional  of  the  Electron  Density.”  Physical 
Review,  vol.  B  45,  p.  13196, 1992. 

Yang,  S.  H.,  M.  J.  Mehl,  and  D.  A.  Papaconstantoupolos.  “Application  of  a  Tight-Binding 
Total-Energy  Method  For  Al,  Ga,  and  In.”  Physical  Review,  vol.  B  57,  p.  R2013, 1998. 


20 


NO.  OF 

COPIES  ORGANIZATION 

2  DEFENSE  TECHNICAL 
INFORMATION  CENTER 
DTIC  DDA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BELVOIR  VA  22060-6218 

1  HQDA 

DAMOFDQ 
D  SCHMIDT 
400  ARMY  PENTAGON 
WASHINGTON  DC  203 10-0460 

1  OSD 

OUSD(A&T)/ODDDR&E(R) 
RJTREW 
THE PENTAGON 
WASHINGTON  DC  20301-7100 

1  DPTY  CG  FOR  RDA 

US  ARMY  MATERIEL  CMD 
AMCRDA 

5001  EISENHOWER  AVE 
ALEXANDRIA  VA  22333-0001 

1  INST  FOR  ADVNCD  TCHNLGY 

THE  UNTV  OF  TEXAS  AT  AUSTIN 
PO  BOX  202797 
AUSTIN  TX  78720-2797 

1  DARPA 
BKASPAR 
3701  N  FAIRFAX  DR 
ARLINGTON  VA  22203-1714 

1  NAVAL  SURFACE  WARFARE  CTR 

CODE  B07  J  PENNELLA 
17320  DAHLGREN  RD 
BLDG  1470  RM  1101 
DAHLGREN  VA  22448-5100 

1  US  MILITARY  ACADEMY 

MATH  SCI  CTR  OF  EXCELLENCE 
DEPT  OF  MATHEMATICAL  SCI 
MADN  MATH 
THAYER  HALL 
WEST  POINT  NY  10996-1786 


NO.  OF 

COPIES  ORGANIZATION 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRLDD 
J  J  ROCCHIO 
2800  POWDER  MILL  RD 
ADELPHI MD  20783-1 197 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  CS  AS  (RECORDS  MGMT) 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 145 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  Cl  LL 
2800  POWDER  MILL  RD 
A  DELPHI  MD  20783-1145 


ABERDEEN  PROVING  GROUND 

4  DIR  USARL 

AMSRL  Cl  LP  (BLDG  305) 


21 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


2  NORTHWESTERN  UNIV 
G  B  OLSON 

MATERIALS  SCIENCE  & 
ENGRG  DEPT 
2225  N  CAMPUS  DR 
EVANSTON  IL  60208 

2  MIT 

DEPT  OF  NUCLEAR  ENGRG 
S  YIP 

77  MASSACHUSETTS  AVE 
CAMBRIDGE  MA  02139 


AMSRLWMMD 
WROY 
AMSRLWMT 
B  BURNS 
AMSRLWMTC 
WDEROSSET 
AMSRLWMTD 
W  BRUCHEY 


ABERDEEN  PROVING  GROUND 


26  DIRUSARL 
AMSRLCI 
C  NEETTBICZ 
AMSRL  Cl  C 
WSTUREK 
AMRSLCIH 
AMSRL  WMB 
A  HORST 
E  SCHMIDT 
AMSRL  WM  BC 
P  PLOSTINS 
AMSRL  WM  BD 
C  F  CHABALOWSKI 
B  FORCH 
BMRICE 
AMSRL  WM  BE 
P  CONROY 
AMSRL  WMM 
D  VIECHNICKI 
GHAGNAUER 
J  MCCAULEY 
AMSRL  WMMA 
RSHUFORD 
AMSRL  WM  MB 
B  FINK 

AMSRL  WMMC 
J  BEATTY 
R  ADLER 
T  HYNES 
J  HIRVONEN 
J  MONTGOMERY 
JDEMAREE 
M  WELLS 


22 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  tor  this  collection  ot  Information  Is  estimated  to  average  1  hour  per  response.  Including  the  time  for  reviewing  Instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  ot  Information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  Including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 

1 .  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

September  1999  Final,  Oct  88  -  Sep  89 

4.  TITLE  AND  SUBTITLE 

A  New  Virial-Theorem-Based  Sevni-Ab-Initio  Method  for  Atomistic  Simulations 

5.  FUNDING  NUMBERS 

L1FX01A 

6.  AUTHOR(S) 

Genrich  L.  Krasko 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory 

ATTN:  AMSRL-WM-MC 

Aberdeen  Proving  Ground,  MD  21005-5069 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ARL-TR-2052 

9.  SPONSORING/MONITORING  AGENCY  NAMES(S)  AND  ADDRESS(ES) 

1 0.SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words ) 

A  new  semi -ab-initio  method  for  atomistic  simulations  based  on  the  virial  theorem  has  been  suggested.  The  method 
is  completely  within  the  realm  of  the  density-functional  theory  and  uses  the  so-called  ‘'reduced"  electron  spin-density 
functional  (SDF).  The  crucial  component  of  the  method  is  the  ansatz  expressing  (for  both  one-component  and 
two-component  systems)  the  electron  density  at  point  r  as  a  superposition  of  "atomic"  densities  due  to  the  neighboring 
atoms.  The  total  energy  of  this  system  is  shown  to  consist  of  three  terms.  The  first  depends  only  on  the  simulation 
volume  and  is  independent  of  the  atomic  configuration.  The  second  and  third,  like  the  embedded-atom  method  (EAM), 
are  the  interatomic  pairwise  interaction  energy,  and  an  ,fN-electron"  term,  which  cannot  be  expressed  as  an  interatomic 
interaction;  it  originates  from  the  electron-correlation  interaction.  The  atomic  densities  are  constructed  using  a  set  of 
polynomial-exponential  functions  resulting  in  an  analytic  form  for  the  pair  interatomic  potential.  The  coefficients  in  the 
atomic  density  expressions  are  found  using  a  calibration  procedure  based  on  performing  a  series  of  ab-initio  calculations 
for  a  few  crystal  modifications  of  this  system.  Success  will  depend  on  whether  the  charge  density  in  a  low-symmetry 
system  under  simulation  will  also  be  close  to  the  true  density  obtainable  from  a  meaningful  ab-initio  calculation;  then, 
the  simulation  results  would  be  identical  to  those  of  the  corresponding  ab-initio  calculation.  To  what  extent  the 
superposition  ansatz  satisfies  this  condition  is  unclear  without  experimental  confirmation. 

14.  SUBJECT  TERMS 

virial  theorem,  atomistic  simulation 

15.  NUMBER  OF  PAGES 

26 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 

OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

UNCLASSIFIED  UNCLASSIFIED  UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

UL 

Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-1 8  298-1 02 


NSN  7540-01-280-5500 


23 


Intentionally  left  blank. 


24 


USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1 .  ARL  Report  Number/ Author  ARL-TR-2052  (Krasko  fPOC:  AdlerV) _ Date  of  Report  September  1999 _ 

2.  Date  Report  Received _ _ _ _ _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for  which  the  report  will 

be  used.) _ _ _____ _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate - 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to  organization, 
technical  content,  format,  etc.) _ _ _ 


Organization 

CURRENT  Name  E-mail  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 

City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the  Old 
or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 


